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Abstract. We show that for large couphng delays the synchronizability of 
delay-coupled networks of identical units relates in a simple way to the spec- 
tral properties of the network topology. The master stability function used 
to determine stability of synchronous solutions has a universal structure in 
the limit of large delay: it is rotationally symmetric around the origin and 
increases monotonically with the radius in the complex plane. We give details 
of the proof of this structure and discuss the resulting universal classification 
of networks with respect to their synchronization properties. We illustrate this 
classification by means of several prototype network topologies. 
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1. Introduction 

Synchronization phenomena of coupled nonlinear oscillators are omnipresent and 
play an important role in physical, chemical and biological systems [H [3] H]. 
Understanding the synchronization mechanisms is crucial for many practical appli- 
cations. 

One of the most intriguing effects occurring in coupled nonlinear systems is the 
synchronization of chaotic dynamics [5^. The notions of chaos and synchronization 
apparently contradict each other. A system is chaotic if small perturbations of 
the system's initial condition are amplified resulting in an unpredictable dynamical 
behavior. Stable synchronization of two systems on the other hand occurs when de- 
viations between the system states decay with time (negative transversal Lyapunov 
exponent). 

The contradiction between these two characteristics is only apparent because the 
decay and amplification occur in different directions in phase space. Perturbations 
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within the synchronization manifold (SM), i.e., the manifold on which the states 
of the two systems are identical, grow due to a positive Lyapunov exponent within 
this manifold giving rise to the chaotic dynamics. On the other hand, perturbations 
transversal to the SM, associated with deviations between the two systems, decay 
due to a negative transverse Lyapunov exponent, thus leading to stable synchroniza- 
tion. Another aspect is the influence of delayed coupling, since time delay renders 
the phase space infinite-dimensional and may on one hand induce instabilities and 
bifurcations, or may on the other hand stabilize unstable states O |7l El |9] . 

Semiconductor lasers are of particular interest in the study of chaos synchroniza- 
tion. The synchronization properties may lead to potential applications, e.g., to new 
secure communication schemes [T0l[IIl[Il[l3[ll[I3[lS[T7l[Tl^ However, it is 
impossible to synchronize chaos in two delay-coupled systems without self-feedback 
for large delays because the synchronized solution is always transversely unstable. 
In coupled lasers this effect leads to spontaneous symmetry breaking, and only 
generalized synchronization of leader-laggard type occurs [20j . 

However, chaos synchronization of two delay coupled systems can be stable if 
each system has self-feedback. For semiconductor lasers this has been realized with 
a passive relay in form of a semitransparent mirror or an active relay in form of 
a third laser in between the two lasers [21\ l22l \23\ [24] . These structures are thus 
interesting for chaos based communication systems. 

For practical applications it is not only necessary that the synchronized solution 
is (linearly) stable, but it is also important how robust the synchronization is to 
noise. Here, nonlinear effects may play an important role. In particular, for syn- 
chronized chaotic systems bubbling [25", "26j plays a key role in this context. This 
effect may lead to occasional desynchronization even for arbitrarily small noise am- 
plitudes. Bubbling has been observed for example in optical [27l |28j and electrical 
[29j systems. This paper reviews and extends our recent letter [30] on synchroniza- 
tion in delay-coupled networks with large delay, providing a universal classification 



2. Structure of the master stability function for large delay 

To determine the stability of a synchronized state in a network of identical units, 
a powerful method has been developed [32, 33j called the master stability function 
(MSF). The MSF for networks with coupling delays has been the subject of recent 
works [331 inSl inSI • Time delay effects play an important role in realistic networks. 
For example, the finite propagation time of light between coupled semiconductor 
lasers [37] [Ml ISl [HI [40l SIl [12 [33 [44^ significantly infiuence the dynamics. Similar 
effects occur in neuronal [33[3a[33l3Hl[39l[5Ql|5ll[52l[53] and biological [54J 
networks. 

In this work we show ^30j that the MSF for networks with coupling delays has a 
very simple structure in the limit of large delays. This allows us to prove a number 
of general statements about the synchronizability of such networks. 

We will first discuss the MSF approach and since we are interested in delay 
coupled systems, we will do this in the context of networks with delay [34 1 [35 l [36], 
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3. Stability of synchronized solutions 

Consider a system of N identical units connected in a network with a coupling 
delay r [35j 

(3.1) j^x\t) = f [x\t)] ^Y.9^jh [x^{t - r)] 

with X* G W^. Here, gij G M is the coupling matrix determining the coupling 
topology and the strength of each link in the network, / is the (non-linear) function 
describing the dynamics of an isolated unit, and is a possibly non-linear coupling 
function. A synchronized solution can only exist if the row sum of the matrix is 
the same for each row, i.e., a = Y^^=i9i3 independent of i. In this case if the 
systems start in a synchronized state = x'^ = - - - = x^ = x(t), they will remain 
synchronized since the feedback term will be equal for all x^ . The synchronized 
solution x{t) is then determined by 

(3.2) ^^^it) = f[^i^t)]+ah[x{t-T)]. 

To calculate the stability of this synchronized solution, we consider the evolution 
of small perturbations ^*(t) 

x\t)=x{t)^C{t). 



Inserting this ansatz into Eq. (3.1) and linearizing in we find 

. N 
(3.3) -Cit) = Df [x{t)] e{t) -^^g^JDh [x{t - r)] ^(t - r), 

where Df and Dh are Jacobians. Using the vector 
the system of linear equations (|3.3|) can be written as 



(3.4) jE{t) = lN0Df [x{t)] E{t) ^g^Dh [x(t - r)] S(t - r), 

where In denotes the A^-dimensional identity matrix and g = {gij) is the coupling 
matrix. We assume that the coupling matrix g is diagonalizable, i.e., there exists a 
unitary transformation U such that 

UgU~^ = diag(cr, 71, 72, ... , 7iv-i)- 

Here, a is the row sum of ^, which is always an eigenvalue of g corresponding to the 
eigenvector (1, 1, . . . , 1). We call this the longitudinal eigenvalue of g. The other 
eigenvalues we then call the transversal eigenvalues of g. 

Diagonalizing g in Eq. (3.4) with the transformation U does not affect the term 
In Df[x{t)] in Eq. (3.4), since this transformation acts only on the matrix In- 
Thus after the diagonalization Eq. (|3.4|) is transformed into N equations 



(3.5) |c(i) = Df [x{t)] m +<yDh [x{t - t)] at - t), 

(3.6) |c(i) = Df [x{t)] m + Ik Dh [x{t - r)] - r) 
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with k = 1,...,A^ — 1. The first equation corresponds to perturbations in the 
direction of the vector (1, 1,..., 1), which act equally on each individual system 
and thus do not cause desynchronization. A growing perturbation in this direction 
indicates that the synchronized solution of the network is chaotic. 

The A^ — 1 other equations in (3.6) describe perturbations transversal to the SM. 
The synchronized solution is stable if and only if these perturbations decay, i.e., if 
the maximum Lyapunov exponent arising from the variational Eq. (3.6) is negative 
for all transversal eigenvalues jk- 

Pecora and Carroll [5j defined a function Amax : C ^ M, which maps a com- 
plex number re'^^ to the maximum Lyapunov exponent arising from the variational 
equation 

= Df [x{t)] m + re"f Dh [3?(t - r)] ^(t - r). 

This function is called the master stability function and it can be calculated nu- 
merically. Once this function is known on a sufficiently large domain in C, one 
can immediately decide for any network structure whether synchronization will be 
stable or not. One only needs to evaluate the MSF at the transversal eigenval- 
ues 7/e of the particular network's coupling matrix. Thus the condition for stable 
synchronization is then 

(3.7) A^ax(7fc)<0, A: = 1,...,7V-1. 

This way the problem has been separated into one part that only depends on the 
dynamics of the individual system, and another part that only depends on the 
coupling topology. 

We will now restrict our analysis to maps, but all ingredients of our arguments 
are also valid for flows and we will point out where the results differ slightly for 
flows. Delay coupled maps have been widely studied because they show similar 
behavior as delay differential equation and interesting synchronization phenomena 
have been found in these systems [55] . 

For delay coupled maps the dynamics in the SM is governed by the equation 

(3.8) Xfe+i = f{xk) + ah{xk-r) 

with r G N and Xk e or eR"^ and the MSF is calculated from 

(3.9) a+i = Df{xk)ik + re'^Dh{xk-r)ik-r- 

Whether the synchronized dynamics is chaotic or not depends on whether the MSF 
evaluated at the eigenvalue re*^ = a, which corresponds to perturbations parallel 
to the SM, is positive or not. 

With the matrix coefficients := Df{xk)^ and := Dh{xk-r) the variational 
equation is given by 

(3.10) a+i = Ak^k + re^^^fc^-r. 

Note that when the delay is changed the dynamics in the SM changes, too. 
Hence, we are not able to make predictions about what happens as r is changed. 
However, at a fixed large value of the delay time r we can compare the Lyapunov 
exponents arising from different values of re*^ in Eq. (3.9). 

We will now analyze the Lyapunov exponents arising from Eq. ( 3.1Q[ ) in the limit 
of large r. We do this in the following steps: first we analyse the two simpler cases, 
where the dynamics in the SM is a fixed point (FP) or a periodic orbit (PO). Then 
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to expand the results to chaotic trajectories x/c in the SM we use the fact that POs 
are dense in a chaotic attractor. 



3.1. A fixed point in the synchronization manifold. For FPs of delay differ- 
ential equations there exists a scaling theory for the FP's eigenvalues in the limit of 
large delay f56| l57l ISH j l59[ l6Q ] . Recently this theory has been generalized to the scal- 
ing of Floquet exponents [61] . In both cases the eigenvalues or Floquet spectrum 
consist of two parts: a strongly unstable part arising from unstable eigenvalues of 
the system without delay and a pseudo-continuous spectrum, for which the real 
part of the eigenvalues approach zero in the limit of large delay. This scaling the- 
ory has been developed for flows. Since we restrict ourselves to maps, we want to 
discuss the main ideas of the scaling theory for maps now. However, each step can 
be done in the same way for flows by applying the large delay theory developed in 

[51 Ell [581 Ei EI]. 

Let us first consider the case, where the dynamics in the SM is a FP, i.e., a 
period T = 1 orbit. In this case the coefficient matrices in Eq. (3.10) are constant 
A = Ak and B = Bk. 

Making the ansatz = ^^^o, we find an equation for the multipliers z 

■re'^^Bz- 



(3.11) 



x{z)=det[A-zI- 



0, 



where / denotes the identity matrix. 

For the strongly unstable spectrum we suppose there is a solution with \z\ > 1. 
Then in the limit of r ^ oo Eq. ( |3.11| ) becomes 

(3.12) det[A-zI] =0. 

Thus in the limit of large delay the eigenvalues z of A with |^| > 1 are also solutions 
of Eq. KU\. 



We are now interested in the pseudo-continuous spectrum, i.e., in the solutions 
with \z\ ~ 1 in the limit of large r. We make the ansatz z = {1 



limit r ^ oo we have (1 + S/r) 
r ^ oo Eq. (3.11) becomes 



e 



(5/r)e^^. In the 
\ and (1 + (5/r) 1. Thus in the limit 



(3.13) = det[A - le'^ + re-^e^^^-^^^] 

with (j) = ujT. The complex equation ( |3.13| ) has three variables 0, a;, and S. Hence, 
it determines real- valued functions 6{uj). As a result, it determines the multipliers 
z{uj). Let us show how S{uj) can be calculated. In particular, we will see that S{uj) 
is independent on the phase ip of the variational equation. 

If B is invert ible, we can calculate the eigenvalues /i 
lowing equation 

= det[-B-\A 



(3.14) 



./e-)-/i]. 



which is a polynomial in /i. This polynomial has exac tly d roots /ij (j = 1, . . . , (i), 
which are eigenvalues of —B~^{A — le^^). Since Eq. (3.14) is independent of ?/^, /i 
does not depend on ^jj as well. 



If B is not invertible, Eq. (3.13) still gives a polynomial in /i, for which the roots 



can be calculated. Then each eigenvalue /i will be a function of uj and one can find 
the branches 



5{uj) = In 



■ In \ fi{oj) \ + Inr. 
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Figure 1. Pseudo continuous spectrum 5{uj) (lines) and location 
of the exact roots (crosses) for the example of a one dimensional 
complex map for r = 3.3 > ro = 3 and r = 2.7 < ro = 3. Parame- 
ters: a 0.4, b 0.2, = 0, r = 30. 



The function /i(cj) can admit the zero value at some point ujq^ i.e., /i(cjo) = 0, in 
the case when the matrix A has an eigenvalue with \z\ = 1. Indeed, as follows from 



Eq. (3.14), for /i = 0, uj = cuo and det B we have 

det[A - le'^'] = det [A - Iz] = 0. 
In all other cases, with det 5 ^ and \z\ ^ 1, the function |/i(co')| is bounded 

< /io < < Ml- 

If there are no strongly unstable eigenvalues, the sign of S determines the stability 
in the limit of large r. It is clear, that 6 increases monotonically with increasing r 
and in particular S is negative for small r and positive for large r. Thus there is a 
minimum radius ro, for which the first eigenvalue branch becomes unstable (5 > 
and thus the MSF changes sign. 

Note that we have obtained the function S{uj) on which the solutions lie in 
the limit of large r but not yet the exact values of u. These are not important 
in the limit of large delay, since the eigenvalues become very dense on the curve 
S{uj). Indeed, the exact values of uj can be calculated from the expression /i(cj) = 
^g-5(u;)gi(^-u;r)^ which implies 

(3.15) Arg /i(cj) = - CJT + 27rk 

for any integer k. Since /i(co') is a known eigenvalue of the matrix —B~^{A — 



/e*^), Eq. (3.15) can be considered as a transcendental equation for determining 



the solutions u = ujk- In particular, Eq. (3.15) implies that the distance between 
the neighboring solutions uJk and uok-i reads 

uJk - ujk-i = - [Arg/i(cjfe-i) - Arg/i(cj/e)] + ^ 
r r 

= 27r/r + O (l/r^) . 

Thus it is proportional to 1/r and the curve 5{uj) is filled densely with roots as 
r ^ oo. 
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Figure 2. Master stability function for a one- dimensional map 
with fixed point dynamics in the synchronization manifold. The 
red regions correspond to A^ax > (synchronized state is unsta- 
ble). The gray regions correspond to Amax < (synchronized state 
is stable). The blue circle indicates the stability boundary given 
by ro according to Eq. (3.16). Already for relatively low values of r 



the blue line matches very well the numerically obtained boundary. 
Parameters of the variational equation: a = 0.4, h = 0.2, r = 20. 




Figure 3. Boundary r(?/^) of stability domain in polar coordi- 
nates for different values of r in logarithmic scale. With increasing 
r r{il)) ro- Parameters of the variational equation: a = 0.4, 
b = 0.2. 



Note that the curve 6{uj) is determined in the bounded interval G [0, 27r] in 
contrast to the case of delay differential equations [57J, where uj varies on the whole 
axis (— oo, oo). 
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The simple case is that of a one dimensional {d = 1) complex map, where 
A = a G C and B = b G C with \a\ < 1. In this case we can explicitly calculate 

S{uj) = \n{\rb\/\a-e''^\). 

For r < (1 — |a|)/|&| all the eigenvalues approach magnitude 1 from the stable side 
and for r > (1 — there are always weakly unstable eigenvalues. Thus the 

MSF changes sign at 

(3.16) ro = (l-|a|)/|6|. 

The pseudo-continuous spectrum for these two cases is depicted in Fig. [T] The 
corresponding MSF is shown in Fig. [2j As r ^ oo, the Amax = contour line 
approaches the circle with radius tq. This is depicted in Fig. [3j where the angle- 
dependence of the critical radius is shown for different values of r in a logarithmic 
scale. For small values of r, the critical radius has a strong angle-dependency. 
However, already for a value of r = 20, the rotation symmetry is almost perfect 
(see Fig.[2|. 

3.2. A periodic orbit in the synchronization manifold. Now consider the 
map 

(3.17) a+i = ^fc a + re'^^Bk ^k-r. 

where A/^ and B}^ are periodic with period T, corresponding to a PO in the SM. 
We consider the case of large delay r ^ T. 

Making a Floquet-like ansatz = Qk^ where qk is T periodic we find 

(3.18) z qk+i = Mqk^ re'^^Bk z~^qk-n 

with n = rmodT G {0, 1, . . . , T - 1}. 

For the strongly unstable spectrum again suppose there is a solution with \z\ > 1, 
then in the limit r ^ oo the term z~^ vanishes and we find 

(3.19) zqk^i = Akqk- 



Using the periodicity of q^^ Eq. (3.19) implies 

T 



det[z^- Afc] =0, 



/e=l 



where z^ is a Floquet multiplier of the system (^k+i = ^fe^fe without delay. 



Hence, if z^ is a Floquet multiplier of Eq. (3.19), with \z\ > 1, then in the limit 



r — > oo there exists also a solution of Eq. (3.17), which is close to z and vice versa. 



The rigorous proof of this fact as well as the convergence for the pseudo-continuous 
spectrum is more tedious and we do not discuss it here. The proof for the case of 
fiows can be found in Ref. [62j. 

For the pseudo-continuous spectrum we again make the ansatz z = {1 -\- S/r)e'^^ 



and taking the limit r ^ oo Eq. (3.18) becomes 

(3.20) e'^ qu+i =Auqk^ re'' e'^^'^^ B^ qk- 
with (j) = ujT. One thus has to solve 

(3.21) = + A + ^B]q = 0, 
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where A = diag{Ai, . . . , At}, 



J = 



/ 
/ . 

/ 



B = 



fji = 7-g-^gH^-0)^ and q = (gi, . . . , qt). The position of the diagonal hnes in the 
matrix B depends on the value of n = rmodT. Taking the determinant of the 
matrix in Eq. ( |3.21 ) results in a polynomial in /i = re^'^e*^^"'^^ (of maximum order 
d X T). Again, the roots /i are functions of u and we can calculate the branches 
S{uj) = — ln|/i(a;)| + Inr, where ip and (j) drop out. As in the case of FPs, one 
can show that the function |/i(cj)| is bounded < jno < |/i(co')| < /ii unless the 
instantaneous system has a Floquet multiplier z with \z\ = 1. 

We have again found the same structure of the MSF : The MSF is rotationally 
symmetric in the complex plane about the origin. If without feedback (r = 0) the 
MSF is positive, then it is constant in the limit of large delay. Otherwise it is a 
monotonically increasing function of r and there is a critical radius ro where it 
changes sign. 



3.3. A chaotic attractor in the synchronization manifold. In every chaotic 
attractor there is an infinite number of unstable periodic orbits (UPOs) embedded. 
It has been long known that the characteristic properties of the chaotic system can 
be described in terms of these PO [63j. Intuitively, the chaotic trajectory follows 
the UPOs closely and "switches" between them, thus averaging over the UPOs in 
the appropriate way allows us to calculate statistical properties of the attractor. 
One of the most important examples is the natural measure of the chaotic attractor, 
which is concentrated at the UPO (hot-spots) and can in fact be expressed in terms 
of the orbit's Floquet multipliers jM } [55 ] , 



Lyapunov exponents arising from variational equations such as Eq. (3.10) have 
been discussed in the framework of periodic orbit theory [66, £71 [63l [68j, too. In 
particular it has been shown [69j that a chaotic attractor in an invariant manifold 
loses its transversal stability in a blow-out bifurcation when the transversely un- 
stable orbits outweigh the transversely stable orbits. To be precise, we divide the 
orbits into two groups of transversely stable and unstable orbits and define [69j the 
transversely stable weight and the unstable weight A^ as 



(3.22) 



T 



where the sum goes over all transversely unstable and transversely stable 
orbits of period T, respectively. Here, jJ^rij) is the weight of the jth orbit, cor- 
responding to the natural measure of a typical trajectory in the neighborhood of 
the jth orbit and XriJ) is the transversal Lyapunov exponent of this jth orbit. 
The weight of a PO is inversely proportional to the product of its unstable Floquet 
multipliers [64^ £8j . The attractor is transversely unstable if and only if in the limit 
of large T 



(3.23) 



A^ > 



A? 



10 



V. FLUNKERT, S. YANCHUK, T. DAHMS, AND E. SCROLL 



Table 1. Stability of chaotic and non-chaotic synchronized so- 
lutions for the three types of networks. 





chaotic dynamics 


PO or FP in the SM 




in the SM (ro < 


(kl<ro) 


(A) iTmaxI < \o\ 


synchr. stable if 


synchr. stable 




|7max| < 




(B) |7max| = kl 


synchr. unstable 


synchr. stable 


(C) iTmaxI > kl 


synchr. unstable 


synchr. stable 






if |7max| < ro 



We now make the connection to the scaling theory for large r. Starting from 
r = (no feedback) transversal Lyapunov exponents XtU) of each orbit can only 
increase with increasing r, as shown above. In particular for large enough r the 
orbits become transversely unstable: either they are already unstable for r = and 
thus remain unstable or the pseudo-continuous spectrum goes to zero and for large 
r it does so from the unstable side. Thus there exists a minimum radius ro, for 
which the condition (3.23) on the weights is fulfilled. Note that since we consider 
the limit r ^ oo we can evaluate Eq. (3.23) at arbitrarily large T, although it is a 
common result of PO theory that formulas such as Eq. (3.22) converge quickly. 

Thus in summary the MSF has the same structure as for FPs and POs (the 
rotation symmetry follows from the rotation symmetry of each Xrij))' 



4. Consequences for synchronization of networks 

Let us now discuss what the structure of the MSF means for the synchroniz- 
ability of networks. We can categorize networks into three types depending on the 
magnitude of the largest transversal eigenvalue 7max in relation to the magnitude 
of the row sum a: (A) the largest transversal eigenvalue is strictly smaller than the 
magnitude of the row sum (|7max| < (B) the largest transversal eigenvalue has 
the same magnitude as the row sum (|7max| = Wl), and (C) the largest transversal 
eigenvalue has a larger magnitude than the row sum (|7max| > Wl). 

At r = |cr| the MSF is positive (ro < |cr|) for chaotic dynamics in the SM and 
negative (|<j| < ro) for dynamics on a stable PO or a FP. This gives us a lower 
or an upper bound on ro and we can thus give the classification as shown in Table 
[l] In networks of type (A) and (B) synchronization on a FP or a PO, which is 
stable within the SM, is always stable. For type (C) this dynamics may be stable 
or not depending on the particular network (value of |7max|) and the dynamics in 
the SM (value of ro). On the other hand chaos synchronization is always unstable 
in networks of type (B) and (C) and it may be stable or not in networks of type 
(A) again depending on the particular network and the dynamics. 

Note that for autonomous flows with a stable PO in the SM we always have 
ro = |cr|, due to the PO's Goldstone mode. Thus for this case we cannot decide 
whether synchronization for type (B) networks will be stable or not. This depends 
on whether the A^ax = contour line of the MSF approaches the circle with radius 
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a) 



(0> 





b) C(Oo 




Figure 4. (a) unidirectionally coupled rings without feedback: al- 
ways class (B) (no chaos synchronization), (b) unidirectionally cou- 
pled rings with feedback (for positive coupling): class (A) (chaos 
synchronization possible). 

ro = |cr|, locally, at the transversal eigenvalues with |7/c| = |cr|, from the outside 
(stable) or from the inside (unstable). 

We now list some examples for the three types of networks (Fig.|4|. The catego- 
rization follows from the eigenvalue structure (spectral radius) for the corresponding 
matrices, which can, for instance, be derived using Gerschgorin's theorem. 

• Mean field coupled systems (all-to-all coupling) have 7/c = for all k and 
are thus of type (A). 

• Networks with only inhibitory connections (negative entries) or only excita- 
tory connections (positive entries) are up to the row sum factor stochastic 
matrices, i.e., the coupling matrix G can be written as 



where P is a stochastic matrix (non- negative entries and unity row sum). 
For stochastic matrices it is well known that the spectral radius is one, 
i.e., all eigenvalues have magnitude smaller than or equal to one. The 
proof utilizes Gerschgorin's theorem [70j. Thus it follows for G that no 
eigenvalues has magnitude larger than |cr| and these networks are of type 



• Any network with zero row sum (<j = 0) is of type (B) (trivial case) or (C). 

• Rings of uni-directionally coupled elements and two bidirectionally coupled 
elements are of type (B) [36j. 

• Unidirectionally coupled rings with self-feedback are for positive couplings 
of type (A). 

• Bidirectionally coupled rings with even number of elements without self- 
feedback are of type (B). 

• Bidirectionally coupled rings with odd number of elements are of type (A). 

In the literature there is a great amount of material on the relation of the spectral 
radius and the row sum for certain types of matrices. These results are immediately 
applicable to our classification. For a concrete network topology the classification 
is of course very simple. 

Networks with a = belong to class (B) (trivial case) or to class (C). This 
confirms the conjecture stated in ^|: Networks for which the trajectory of an 
uncoupled unit is also a solution of the network (a = 0) cannot exhibit chaos 
synchronization for large coupling delay. 



G = aP , 



(A) or (B). 
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Figure 5. Master stability function (shown in color) for delay 
coupled logistic maps (Eq. (4.2)). Panel (a) corresponds to A = 3.2 
and the stable period-2 orbit within the SM. Panel (b) corresponds 
to A = 3.8 and a chaotic attractor within the SM. In both cases 
the delay is chosen as r = 30. 



For the chaotic case there may exist another radius r5, with < < ro, where 
the first PO in the attractor loses its transverse stability and the attractor undergoes 
a bubbling bifurcation [26l |7T| |42] . Then any network with < |7max| < 
will exhibit bubbling in the presence of noise (or parameter mismatch), while any 
network with |7max| < ^6 will show stable synchronization, even in the presence of 
noise. 

In order to illustrate the obtained results, let us consider the following example 
of linearly coupled logistic maps 

N 

(4.1) X^^, = XX^{1 - X^) ^Y.3mjXi_^ 

J = l 

with the zero row sums a = 0. The MSF is calculated from the following delayed 
system 



(4.2) a+i = A(l-2xfc)a + re^^e 



k — T") 



where the dynamics on the synchronization manifold x/c is determined by the map 
Xk+i — Ax/e(l — x/e). Figure [s] shows numerically computed MSF, i.e., the largest 



Lyapunov exponent of the system (4.2) for two different cases: A = 3.2 and A = 3.8, 
which correspond to the stable period-2 state and chaos, respectively. The delay 
is set to r = 30. In both cases, the MSF are radially symmetric. In the stable 
periodic case (panel (a)), there exists a critical radius ro where the MSF changes 
sign, which determines the synchronizability of a given coupled system. In the 
chaotic case (panel (b)) the MSF is close to a positive constant, i.e., any coupling 
configuration will be unstable. 

For the parameters used in panel (b) the dynamics exhibits so-called strong chaos 
[72j . Strong chaos is characterized by a Lyapunov exponent that stays constant with 
increasing delay time. As recently shown [72j chaos synchronization is not possible 
in the large delay limit for strong chaos. The other case of weak chaos occurs 
when the largest Lyapunov exponent scales as 1/r for r ^ oo. In this case chaos 
synchronization is possible and the critical radius is determined by 

(4.3) ro = \a\e-^-^ , 
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where is the maximum Lyapunov exponent of the system [72\. 

5. Experimental setup for finding the critical radius 

We now propose an experimental method for determining the critical radius tq. 
Consider two elements coupled in the following network motif 

4+1 = + ^K^l-r) + ^H^l-r)^ 

4+1 = /(4) + ^H^I-t) + ^K^k-r)^ 

where ji and v are the self feedback strengths and the coupling strengths, respec- 
tively. Suppose we are able to change the self-feedback strengths ji and the coupling 
strengths for example by using gray filters in an optical experiment. 
Let us choose 

/i=^(cr + r) and h'=\{cF — r). 
Then the dynamics in the SM is given by 

Xk+i = f{xk) + crh{xk-r), 
while the variational equation transverse to the SM is given by 

^/e+l = Df{Xk)^k + rDh{Xk-r)ik-T' 

Thus by changing r (for fixed a) and checking whether the two elements synchronize 
we are able to probe the MSF along the real axis at the radius r. Due to the 
monotonicity we can use a root-finding algorithm such as the bisection method 
to find ro to high accuracy with little iterations of the experiment and without 
knowledge of the functions / and h. We can repeat this procedure for other values of 
<j and obtain the critical radius as a function of ro (cr) . Thus from this rather simple 
setup we can decide for any network of these elements whether synchronization is 
stable or not. 

As an example we consider two optoelectronically coupled lasers 
|pi=nipx, 

(5.1a) -^^^1 = P + l~^pi(t - r) + vp2(t - r) - ni - (1 + ni)pi, 

d 

— P2 n2p2, 

(5.1b) -^^^2 = P + IJip2{t - r) + vpiit - r) - 712 - (1 + ^2)^2, 

where pi and is the intensity and the carrier density of the ith laser, respectively. 
The pump current of each laser is modulated by the delayed intensities according 
to the coupling scheme depicted in Fig. |6] Such feedback can be realized by using 
photodiodes to measure the intensities of the arriving signals and modulating the 
pump current accordingly. The bidirectional coupling has strength v and the self- 
feedback of each laser has strength p. For this setup we have numerically calculated 
ro(cr) in the same manner as it would be done in an experiment: We choose a value 
of (T, an interval 1^ — [^min? ^max 

1 for the r-domain and an initial value of r = ro- 
We consider the systems to be synchronized, if the relative synchronization error 

(|Pi -P2I) 

' \[{Pl)^{P2)] 
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If 

Figure 6. Schematic setup for determining the critical radius in 
an experiment. 




Figure 7. Numericahy calculated critical radius ro as a function 
of (T (solid curve) for the system of optoelectronically coupled lasers 
corresponding to Eqs. (5.1 ). A network can only have a stable syn- 
chronized solution if the magnitudes of its transversal eigenvalues 
are below the curve. The curve is calculated up to an absolute 
error of 10~^. The dashed line shows the diagonal line ro = cr. 
Parameters: Sq lO"'^, p = 1, T = 200, r 2000. 



is smaller than a threshold Sq. We then simulate the system and use the bisection 
method to find the synchronization threshold ro (up to a desired accuracy) in the 
interval We can then use the calculated value of ro as an initial guess for 
neighboring a- values and thus follow the curve ro(cr). 

The result is depicted in Fig.[7| where the solid curve shows ro(cr) and the dashed 
line corresponds to ro = cr. For small values of a, i. e., weak feedback, the dynamics 
is a PO and due to the Goldstone mode we have ro ^ cr. For larger values of cr, 
the system becomes chaotic and ro < cr. For a given value of cr, a network has 
a stable synchronized solution if and only if all transversal eigenvalues 7/c of the 
corresponding coupling matrix have magnitude |7/c| < Tq{(j). 

A similar method has been implemented in |73j. Here the structure of the MSF 
was confirmed experimentally using optoelectronic oscillators. 

6. Conclusion and outlook 

In conclusion we have shown that the MSF has a simple structure in the limit of 
large delay: it is rotationally symmetric around the origin and either positive and 
constant (if it is positive at the origin) or monotonically increasing and becomes 
positive at a minimum radius ro- This structure allowed us to prove a recent con- 
jecture [35] about synchronizability of chaotic elements. Furthermore, we classified 
network structures into three types depending on the magnitude of the maximum 
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transversal eigenvalue in relation to the magnitude of the row sum and showed 
that these network types have distinct synchronization properties. Using several 
prototype networks like all-to-all or ring topologies we illustrated the scope of these 
three classes. By means of a motif of two coupled nodes with feedback loops, we 
proposed a very simple scheme with which the critical radius Tq can be found ex- 
perimentally. Instead of mapping out the entire domain of the master stability 
function, the knowledge of the rotational symmetry allows the two- node motif to 
predict stability of any network using the same local dynamics in the case of large 
delay times. 

The rotational symmetry of the MSF has previously been found numerically 
[361 152] , In Ref. [52j the same structure of the MSF has been found for a PO in 
the SM for which the period T is approximately equal to the delay time r. For this 
case the structure of the MSF has also been derived analytically in [52j. Note that 
this case is complementary to the situation T <C r that we looked at in this paper. 
So the structure of the MSF that we found seems to be valid in even more general 
cases. 

VF acknowledges financial support from the German Academic Exchange Service 
(DAAD). This work was performed in the framework of SFB 910. 
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